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I. INTRODUCTION 



The ability to compute the details of the gravitational radiation produced by compact astrophysical sources, such as 
coalescing black holes, is of major importance to the success of gravitational wave astronomy. The simulation of such 
' systems by the numerical evolution of solutions to Einstein's equations requires a valid treatment of the outer bound- 
I ary. This involves the mathematical issue of proper boundary conditions which ensure a well-posed initial-boundary 
> value problem (IBVP), the computational issues of the consistency, accuracy and stability of the finite difference 
' approximation, and the physical issue of prescribing boundary data which is free of spurious incoming radiation. The 
1 physical issue can be solved by locating the outer boundary at future null infinity 1+ in a conformally compactified 

■ spacetime . Since no light rays can enter the spacetime through X+ , no boundary data are needed to evolve the inte- 
, rior spacetime. In addition, the waveform and polarization of the outgoing radiation can be unambiguously calculated 
' at 1+ in terms of the Bondi news function. This global approach has been applied to single black hole spacetimes 

Q I using Cauchy evolution on hyperboloidal time slices (see 0, Q for reviews) and using characteristic evolution on null 
O^" hypersurfaces(see Q for a review). However, most current computational work on the binary black hole problem 
jj^ involves Cauchy evolution in a domain whose outer boundary is a finite timelike worldtube. This approach could be 
bi). justified physically if the boundary data were supplied by matching to an exterior solution extending to X+. The 

■ present work is part of a Cauchy-characteristic matching (CCM) project |^|^ to achieve this by matching the interior 
• i-H , Cauchy evolution to an exterior characteristic evolution. The success of CCM depends upon the proper mathematical 

■ and computational treatment of the Cauchy IBVP. A key issue, on which we focus in this paper, is the accurate 
^ , preservation of the constraints in the treatment of the Cauchy boundary. The results have important bearing on the 
- - < matching problem. 

Early computational work in general relativity focused on the initial value problem in which data on a spacelike 
hypersurface S determines the evolution in its domain of dependence. However, in the simulation of an isolated system, 
such as a neutron star or black hole, S typically has an outer timelike boundary B, coincident with the boundary of 
the computational grid. The initial-boundary value problem addresses the extension of the evolution into the domain 
of dependence oi S U B. The IBVP for Einstein's equations is still not well understood due to complications arising 
from the constraint equations. See for a review of the mathematical aspects. Well-posedness of the mathematical 
problem, which guarantees the existence and uniqueness of a solution with continuous dependence on the Cauchy 
and boundary data, is a necessary condition for computational success. A well-posed IBVP for a linearization of the 
Einstein equations was first presented by Stewart , and the first well-posed version for the full nonlinear theory was 
established by Friedrich and Nagy |^ . These works have influenced numerous inv estig ations of the implementation 
of boundary conditions in numerical relativity Hfl IH IH [H HI l^l . However, 

the combined mathematical and computational aspects of the IBVP still present a major impasse in carrying out the 
long term simulations necessary to compute useful waveforms from the inspiral and merger of black holes. 

The IBVP for general relativity takes on one of its simplest forms in a harmonic gauge, in which the well-posedness 
of the Cauchy problem was first established (26l | . In previous work, we formulated a well-posed constraint preserving 
version of the harmonic IBVP and implemented it as a numerical code 6], the Abigel code. The IBVP was demon- 
strated to be well-posed for homogeneous boundary data (or for small boundary data in the sense of linearization 
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off homogeneous data) and the code was shown to be stable and convergent. Nevertheless, numerical evolution was 
limited by the excitation of exponential instabilities of an analytic origin. Insight into the nature of these unstable 
modes and mathematical and computational techniques for dealing with them have been developed in a study of the 
evolution problem in the absence of boundaries. These studies were carried out on toroidal spatial manifolds, equiva- 
lent to the imposition of periodic boundary conditions. In this paper, we apply the techniques which were successful 
in the periodic case to the general harmonic IBVP. Although our applications here are limited to test problems, we 
expect the techniques will be of benefit in furthering recent progress in the simulation of black holes by harmonic 
evolution 25, 27, .22,]. 

Our previous work on the IBVP Q was based upon a reduced form of the harmonic system due to Fock [2^ . We 
present our results here in the more standard form 0, 113 used in most analytic work. We have described and tested 
a finite difference evolution code based upon this standard harmonic formulation in [sij . In Sec. ^ we summarize 
the main analytic and computational features of this evolution code. The formalism includes harmonic gauge forcing 
terms and constraint adjustments. In principle, gauge forcing terms [3^ allow the simulation of any nonsingular 
spacetime region. Gauge forcing not only allow the flexibility to "steer" around pathologies that might otherwise 
arise in standard harmonic coordinates but it also allows universal adaptability to carry out standardized tests for 
code performance, such as the AppleswithApples (AwA) tests The formalism also includes constraint adjustments 
which modify the nonlinear terms in the standard harmonic equations by mixing in the harmonic constraints. Such 
constraint adjustments have proved to be important in harmonic evolution in suppressing instabilities in a shifted 
version of the AwA gauge wave test 31] and in simulating black holes 27, 28i]. 

Constraint adjustments, when combined with a flux conservative form of the equations, can lead to conserved 
quantities which suppress exponentially growing error modes. As a dramatic example, the AwA gauge wave tests for 
the Abigel harmonic code show an increase in error of more than 12 orders of magnitude, after 100 grid crossing times, 
if flux conservation is not employed |l7j ] . The error arises from exponentially growing solutions to the reduced equations 
which have long wavelength and satisfy the constraints, so that it is unaffected by either numerical dissipation or 
constraint adjustment. The underlying conservation laws only apply to the prinicple part of the system. The are 
effective when the nonlinear terms corresponding to the unstable modes are small, or can be adjusted to be small by 
mixing in the constraints. This is the case with the AwA gauge wave. However, as we discuss in Sec. ^ a shifted 
version of the gauge wave test excites a different type of exponential instability which does not satisfy the constraints 
and whose suppression requires constraint adjustment in addition to flux conservation. This raises the caveat that 
the identification of the nonlinear instabilities in the system is critical to the effectiveness of constraint adjustment 
and conservative techniques. 

In See's IIIl| and IIVI we describe how our earlier proof ^ of the well-posedness of the harmonic IBVP for Einstein's 
equations extends to the generalized formulation considered here. Constraint preservation imposes requirements 
coupling the intrinsic metric and extrinsic curvature of the boundary, analogous to the momentum constraint on 
Cauchy data. However, the boundary data has fewer degrees of freedom than the corresponding Cauchy data, as for 
a scalar field where either the field or its normal derivative, but not both, can be specified on the boundary. This 
couples the evolution to the specification of constraint-preserving boundary data in a way that the standard theorems 
used to establish well-posedness only apply to boundary data linearized off homogeneous boundary data. In the case 
of large boundary data, one is only assured that the solution, if it exists, does satisfy the constraints. 

In Sec. we present the techniques used to implement the formalism as a finite difference evolution-boundary 
algorithm. Since the preliminary testing of the Abigel code, considerable improvement has been made in the numer- 
ical techniques. The long term performance of the evolution algorithm in the AwA gauge wave test with periodic 
boundaries has been improved by use of semi-discrete conservation laws 113. The study of a model nonlinear scalar 
wave 34] shows how these semi-discrete conservation laws can be used to formulate stable algorithms for more general 
boundary conditions. 

In Sec. I VII we demonstrate the stability and second order convergence of the code for periodic, Neumann, Dirichlet 
and Sommerfeld boundary conditions. These tests are performed using gauge wave and shifted gauge wave metrics 
for which the exact solutions provide the correct boundary data. In principle, the knowledge of the exact boundary 
data avoids the need for constraint-preserving boundary conditions. However, numerical noise can generate constraint 
violating error. This is investigated by comparing these constraint free boundary algorithms against the constraint- 
preserving algorithm discussed in See's ITTll and lIVl In order to eliminate the complication of sharp boundary points, 
the tests are carried out by opening up the 3-torus (periodic boundary) into a 2-torus times a line, with two smooth 
2-toroidal boundaries. 

Just as the "3-1-1" decomposition — (t, x*) is useful in describing the geometry of a spacelike Cauchy hypersurface 
5 at i = 0, the decomposition — (a;°,a;), with = (t,y,z) is useful in describing the geometry of a timelike 
boundary S at a; = 0. We distinguish between these decompositions by using indices i,j,k, ... near the middle of the 
alphabet or indices a,b,c,... at the beginning of the alphabet. In order to avoid excessive notation we indicate the 
corresponding intrinsic 3- metrics by hij and hat, with inverses /i'^ and Where this might cause confusion, we 
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introduce subscripts, e.g. hs = det{hij) and /ig = det{hab)- Curvature conventions follow |30| . This introduces some 
sign changes from the treatment in [g which was based upon the conventions in 29] . We use the shorthand notation 
daf = f,a where confusion does not arise. 



II. THE GENERALIZED HARMONIC CAUCHY PROBLEM 



We summarize here the treatment of the generalized harmonic Cauchy problem on which the evolution code de- 
scribed in jSlj is based. Generalized harmonic coordinates x" = (t, x*) — {t,x,y,z) are functionally independent 
solutions of the curved space scalar wave equation, 



(2.1) 



where the gauge source terms t^{x°' , gap) can have functional dependence on the coordinates and the metric |32|. In 

(2.2) 



terms of the connection T'^ these harmonic conditions take the form 



where 



C := -r^ = 0. 



and 7''^ = \/-~gg^'^ ■ The standard harmonic reduction of the Einstein tensor (see e.g. 113,123) 



(2.3) 



(2.4) 



where is treated formally as a vector in constructing the "covariant" derivatives V^F". Explicitly in terms of the 
metric. 



ap„/3.pp^p.^ - V~g{dc.g"^)d0g^''' + -^g^P{dpg)dag^'' 



2^gE^^ - da{g''^dp^^'')~2^~gg 



(2.5) 



Thus the standard harmonic evolution equations are quasilinear wave equations E'^'^ = for the components of the 
densitized metric 7^"^. The principle part of l|2.5|l has been expressed in the flux conservative form 



(2.6) 



In the particular case of the AwA gauge wave, the vanishing of the nonlinear source terms then leads to the exact 
conservation laws for the quantitities 



IMS' 



g^'^dar^'dV. 



(2.7) 



The semi-discrete version of these conservation laws suppresses instabilities in the AwA gauge wave test. (See [sTl for 
details.) 

This standard treatment of the Cauchy problem in harmonic coordinates generalizes in a straightforward way to 
include harmonic gauge source terms and constraint adjustments of the form 



(2.8) 



The reduced equations then become 



(2.9) 
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Since the gauge source terms and constraint adjustments do not enter the principle part, (|2.9(l also constitute a system 
of quasilinear wave equations with a well-posed Cauchy problem. 

The solutions of the generalized harmonic evolution system 1)2.9(1 are solutions of the Einstein equations provided 
the harmonic constraints 

r'' - = (2.10) 

are satisfied. The Bianchi identities, applied to H2.9|l . imply that these constraints obey the homogeneous wave 
equations 

V"VaC^ + RfiC - 2V^{CPA^;'') = 0, (2.11) 
where via (|2.9|l the Ricci tensor reduces to 

^ V^^C^) - CP{A^;'' - ^gf^'g^pAf). (2.12) 

The well-posedness of the Cauchy problem for 1(2. 11(1 implies the unique solution = in the domain of dependence 
of the initial Cauchy hypersurface S provided the Cauchy data 'y^'^ls and dt'j'^'^ls satisfy C^ls = dtC^ls = via 
(12.9(1 . It is straightforward to verify that Cauchy data on S which satisfy the Hamiltonian and momentum constraints 
= and the initial condition — also satisfy 9tC'' = on 5 by virtue of the reduced harmonic equations ((2.9(1 . 
In addition to the standard Cauchy data for the 3 + 1 decomposition, i.e. the intrinsic metric and extrinsic curvature 
of S subject to the Hamiltonian and momentum constraints, the only other free data are the initial choices of 7*", 
equivalent to the initial choices of lapse and shift. 

The standard harmonic reduction ((2.4(1 differs from Fock's p9| harmonic formulation adopted in Jj] by the adjust- 
ment 

= ^C<^^g'^)"d^g - ^g^'^rO^g. (2.13) 

In the tests of the boundary algorithms in Sec. IVII we have examined the effects of the adjustments 

= -^C^d^iV^g^^n (2.14) 
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A^" = "'^I^^V ^C^ (2.15) 
= ^^^C^f^V'h, (2.16) 

v-5" 

where the a.; > are positive adjustable parameters and 

epa = 9pa - -^{'^pt)^at (2.17) 

is the natural metric of signature (-1- -I- -I--I-) associated with the Cauchy slicing. The adjustments ((2.14(1 and ((2.15(1 
were effective in suppressing instabilities in the shifted gauge wave tests without boundaries |3lj |. (For numerical 
purposes, the limit as ^ in (12.15(1 can be reg ularized by a small positive addition to the denominator.) The 
adjustment (12.16(1 leads to constraint damping 35] in the linear regime and has been used effectively in black hole 
simulations 28] but was not effective in the nonlinear shifted gauge wave test. 



III. WELL-POSEDNESS OF THE HOMOGENEOUS INITIAL-BOUNDARY VALUE PROBLEM 

In prior work we showed how the well-posedness of the Cauchy problem for the standard harmonic formulation of 
Einstein's equations extends to the IBVP with homogeneous boundary data. Here we show how the well-posedness of 
the IBVP extends to include harmonic gauge forcing terms and constraint adjustments. Our approach is based upon 
a theorem of Secchi 41] for first order, quasilinear, symmetric hyperbolic systems. For that purpose, we recast the 
reduced system ((2.9(1 in first order form in terms of the 50-dimensional column vector u — -^(7"'^, T"^, X"^ , y"^ , Z°'^) 
(the transpose of the row vector ^u) where T"'^ = 9*7"^, A'"'' = 9a;7"'^, 3^"^ = c'y7"'^ and Z"'' = 9^7"''. Equation 
((2.9(1 then has the quasilinear, symmetric hyperbolic form 



A*(9tu -I- A^^iU = S 



(3.1) 
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where A* and A* are 50-dimensional symmetric matrices, with A* positive definite, and S is a 50-dimensional column 
vector, all depending only on the components of u. For explicit values of the matrices A* and A' see App. ^ 

Consider the IB VP for the reduced Einstein equations for the generalized harmonic system in the domain t > 0, 
X < 0, with timelike boundary ;S at a; = 0. (Any boundary can be constructed by patching together such pieces.) 
Secchi's theorem establishes well-posedness of the IBVP for (|3.1I) under weak regularity assumptions by imposing a 
homogeneous boundary condition based upon the energy norm 

S = - / uudxdydz. (3.2) 

and the associated local energy flux across the boundary 

= i'^uA^u. (3.3) 

The chief requirement for a well-posed IBVP is that the boundary condition be maximally dissipative, which allows 
energy estimates to be established. Explicitly, the requirements at a; = are 

• that the kernel of the boundary matrix A^ must have constant dimension, 

• that the boundary condition must take the matrix form Mu = for homogeneous boundary data, where M is 
independent of u and has maximal dimension consistent with the system of equations 

• and, for all u satisfying the boundary condition, that the boundary flux satisfy the inequality 

> 0. (3.4) 



In addition, compatibility conditions between the Cauchy and boundary data at t = x = must be satisfied. 

As a result of reducing the harmonic Einstein equations to first order symmetric hyperbolic form, there are variables 
associated with characteristics which propagate tangential to the boundary. Secchi's theorem applies to this case 
of a characteristic boundary in which the matrix A^ is degenerate. The matrix M can depend explicitly on the 
coordinates but not upon the evolution variables u. The maximal dimension of M is directly related to how many 
variables propagate toward the boundary from the exterior. Variables in the kernel of A^, which propagate tangential 
to the boundary, and variables which propagate from the interior toward the boundary require no boundary condition. 
The possible choices of matrix M correspond to the nature of the boundary condition, e.g. Dirichlet, Sommerfeld or 
Neumann. 

We now formulate such maximally dissipative boundary conditions for the reduced harmonic system. The principal 
part consists of the term g°'^dadf3j'^'''^ , which represents a wave operator acting on the 10 components 7'*". As a 
result, the energy flux equals the sum of 10 individual fluxes, each of which corresponds to the standard energy 
flux constructed from the stress-energy tensor obtained by treating each component of 7'^" as a massless scalar field. 
Thus we can pattern our analysis of the IBVP on the scalar wave equation 

g^'-a^a,* = 0, (3.5) 

which is discussed in App.^ When <I> is regarded as an independent scalar field rather than as a component of 7^"^, 
the linearity of (|3.5|l implies well-posedness for any of the dissipative boundary conditions discussed in App. ^ e.g. 
generalized Dirichlet, Neumann or Sommerfeld boundary conditions. 

We can decompose the boundary conditions for the full gravitational problem into the five 10-dimensional subspaces 
corresponding to 7"^, T"^, X"^ , and Z'^^ . The subspace corresponding to 7"^ lies in the kernel of A^ and 
requires no boundary condition, i.e the boundary values of 7"^^ are determined by boundary values of the other 
variables. The remaining four subspaces require boundary conditions analogous to the scalar field case. The extra 
complications in the quasi-linear gravitational case are to check that the constraints are satisfied and that the linear 
subspace Mu = is specified independent of u. 

For that purpose, we adapt the harmonic coordinate system to the boundary. In doing so, we use the 3-1-1 
decomposition natural to the boundary at x = and write x" = {x°',x), where — {t,y,z). Our approach is 
motivated by the observation that the initial value problem is well-posed so that refiection symmetric Cauchy data in 
the neighborhood —e<x<e produces a solution with even parity. From the point of view of the IBVP in the region 
X < 0, this solution induces homogeneous boundary data of either the Neumann or Dirichlet type on each metric 
component, depending upon how that component transforms under reflection. 
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A. The boundary gauge 

Local to the boundary B, it is always possible to choose Gaussian normal coordinates so that the boundary 

is given by i = and = vanishes at the boundary. Now consider a coordinate transformation to generalized 
harmonic coordinates (x°,a;) satisfying Ox" = — F". Since 17^" = 0, the analysis in Add. IXI shows that it is possible 
to solve this wave equation with either the homogeneous Dirichlet condition $ = or the homogeneous Neumann 
condition — 0. For the harmonic coordinate x, we choose the Dirichlet condition x\b = 0, so that B remains 
located at a; = 0; and for the remaining harmonic coordinates, we choose the Neumann condition dxX°'\B = 0. 

These choices exhaust the boundary freedom in constructing generalized harmonic coordinates. After transforming 
to these harmonic coordinates, the metric satisfies 

7"1b = 0. (3.6) 
We enforce this as the dissipative Dirichlet boundary condition 

r™|B = 0, (3.7) 

subject to initial data satisfying H3.6|l . 

In this boundary gauge, the results of Add. IXIimolv that the homogeneous Neumann boundary conditions 

X^B - = 0, (3.8) 

are also dissipative, as well as the Sommerfeld-like conditions 

{T^^ + X^^)\b^{T''^ + X'''')\b = 0. (3.9) 

Given that the boundary gauge condition (|3.6(l is satisfied, any combination of Dirichlet, Sommerfeld-like and 
Neumann boundary conditions on the remaining components leads to a well-posed IBVP for the reduced generalized 
harmonic system (|2.9f) . Note, however, that a Sommerfeld condition corresponding to the null direction defined by 
normal to the boundary. 



V 9xx Y 9xx 



does not satisfy the technical condition in Sccchi's theorem that M be independent of u. 



B. The homogeneous constraint-preserving boundary condition 

The constraints = F^ — F^ satisfy the homogeneous wave equation (|2.11() , which can be reduced to symmetric 
hyperbolic form. Consequently, we can force them to vanish by choosing boundary conditions for the reduced system 
which force to satisfy a maximally dissipative homogeneous boundary condition. In doing so, we shall require that 
the normal components of the gauge source function and constraint adjustment vanish on the boundary, i.e. 

r^B = 

= 0. (3.11) 

These requirements are a consequence of the reflection properties of the boundary in the homogeneous case. 

We proceed by imposing boundary conditions on the evolution system consisting of the Dirichlet boundary-gauge 
conditions 13.6|l and the Neumann conditions 1)3.8(1 . which imply 

= 9.7""|b = 9.7"% = 9..g|6 = 0. (3.12) 

As we shall show, ((3.6|l and (|3.12(l further imply that the constraints obey the maximally dissipative homogeneous 
boundary conditions 



atC^lB = 
dxClB = 0. 



(3.13) 
(3.14) 
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Note that H3.13|l and H3.14|l are satisfied by any geometry which has local reflection symmetry with respect to the 
boundary if the gauge source functions share this symmetry (as a vector). In order to verify H3.13|l . we note that 

dtC^ = -dt (-^{d^i^^ + dain + r") (3.15) 



vanishes on the boundary as a consequence of (|3.11|) and (|3.12f) . 
Similarly, in order to verify H3.14|l . we note that 



d^r = -d, ^(5,7"" + 967"") + r-^ (3.16) 



reduces on the boundary to 



d,nB = -[^-^dh-- + d.,T-yB- (3.17) 

Now we use the i?"^ — component of the evolution equation to eliminate 3^7"'^. Substitution of the assumed 
boundary conditions into H2.5() gives 

'^"^1- - ^""'^ 527^^18. (3.18) 



Then (jTHJl gives 
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which, together with H3.11|l and H3.17|l . establishes H3.14|l . 

In summary, the maximally dissipative boundary conditions (|3.t)|) and H3.8|) for the evolution variables u imply the 
maximally dissipative homogeneous boundary conditions H3.13|l and (|3.14|) for the constraints. This establishes the 
boundary data necessary to show that the harmonic constraints — Q propagate if the initial Cauchy data satisfy 
the Hamiltonian and momentum constraints. 

The remaining ingredient necessary for a well-posed initial-boundary problem is the consistency between the initial 
Cauchy data and the boundary data. The order of consistency determines the order of differentiability of the solution. 
The simplest case to analyze is prescription of initial data with locally smooth reflection symmetry at the boundary, 
e.g Cauchy data 7"''(0, a;, y, z) for x < whose extension to a; > by reflection would yield a C°° function in the 
neighborhood of a; = 0. Such data are automatically consistent with the homogeneous boundary conditions given 
above. 

We could also have established the boundary conditions (|3.14|) for the constraints indirectly, but with a more 
geometric argument, by noting that the boundary conditions (|3.12() for the metric imply that the extrinsic curvature 
of the boundary vanishes, which in turn implies that the G"*^ components of the Einstein tensor vanish at the boundary 
(see App. IB)| . Thus, along with the evolution equation E"-^ = 0, H2.9|) implies 

- V("C^) -I- ig'^^VaC" + A""" = 0. (3.20) 

This, combined with (|3.11|) - (|3.13|) . now gives H3.14|l . We take this geometrical approach in the next section in 
analyzing the difficulties underlying extending these results to include inhomogeneous boundary data. 



IV. INHOMOGENEOUS BOUNDARY DATA 

In the preceding section we showed that the generalized harmonic evolution system E^^ = gives rise to a well- 
posed constraint-preserving IB VP with a variety of maximally dissipative homogeneous boundary conditions. The 
well-posedness of this IBVP for the reduced system extends by standard arguments to include inhomogeneous data. 
However, only for special choices of the boundary data is the solution guaranteed to satisfy the constraints. For 
example, in the simulation of a known analytic solution to Einstein's equations, when the initial data and boundary 
data are supplied by their analytic values, the well-posedness of the IBVP guarantees a unique and therefore constraint- 
preserving solution. This example is important in code development for carrying out tests of the type reported in 
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Sec. I VII Similarly, when an exterior solution of the Einstein equations is supplied in another patch overlapping 
the boundary, the induced boundary data are constraint-preserving. This possibility arises in Cauchy-characteristic 
matching. Except for such examples, the specification of free boundary data for the reduced system does not in 
general lead to a solution of the Einstein equations. 

We now formulate conditions on the inhomogcncous boundary data which are sufficient to guarantee that the IB VP 
for the generalized harmonic system leads to a solution that satisfies the constraints. However, because the matrix 
M used in formulating these constraint-preserving boundary conditions depends on the evolution variables u, the 
well-posedness of the resulting IBVP follows from Secchi's theorems only for "small" boundary data, in the sense of 
linearization about a solution with homogeneous data. Well-posedness for "large" boundary data remains an open 
question. 

As in the homogeneous case, we consider the IBVP in the domain i > 0, a: < 0. We first generalize the boundary 
gauge condition (|3.()|l . Given harmonic coordinates x" such that ^'^"'\b = at x = 0, consider the effect of a harmonic 
coordinate transformation x^{x'^) satisfying Dx^ — — P'^. For simplicity, we assign initial Cauchy data = x^ and 
dtx^ = dtx^ at t = and require that the boundary remain located at i = 0. By uniqueness of solutions to the IBVP 
for the wave equation, this Dirichlet boundary data for x determines that x = x everywhere in the evolution domain 
so that i'^ = {x°',x). The coordinates are uniquely determined by the Neumann boundary data = dxX°'\B, 
which can be assigned freely. Under this transformation, 7^° = q'^^^^ at the boundary, so that ■y^'^ /-f^^ can be freely 
specified as boundary data (Consistency conditions would require that q"' — at t = x — 0. However, the 

initial value of g"^ can be freed up by a deformation of the initial Cauchy hypersurface in the neighborhood of the 
boundary.) In the remainder of this section, we assume that this transformation has been made and that the 
coordinates have been relabeled a;^. Thus, in this coordinate system, 

xa 

i^i^) = (4-1) 

is free boundary data representing the choice of boundary gauge aX x — Q. 

Although it might seem that such a gauge transformation would not disturb the well-posed nature of the IBVP 
problem for the reduced system, already (|4.1() introduces the complication that Dirichlet boundary data for 7^" 
cannot be determined from (7° unless 7^^ is known on the boundary. But the boundary condition H3.12f) for 7^^ is 
of Neumann type so that 7^^|e can only be determined after carrying out the evolution and, in particular, depends 
on the initial Cauchy data. This is the first source of the u-dependence of the matrix M expressing the boundary 
condition. It would seem that this problem might be avoided by using 7^''/7^^ as an evolution variable but, as will 
become apparent, more complicated (differential) problems of this type arise in introducing inhomogeneous constraint 
preserving boundary data for the remaining variables. 

The complete set of free boundary data for the generalized harmonic system, which in the inhomogeneous case 
replaces the maximally dissipative homogeneous boundary conditions H3.6|l and H3.12|) . consist of 



where, in accord with H4.1|l . 



q- (9^9"^9'^'') - ('7^9'^5^7"^g^5^7'"')b, (4.2) 



q^ = (g", 1) = —Wx = ^=n^ (4.3) 



with the unit outward normal to the boundary. 

Note that a dissipative Neumann condition must be formulated in terms of the normal derivative and not the 
dx derivative used in the previous section where g" — 0. For this purpose, it is useful to introduce the projection 
tensor into the tangent space of the boundary, 

= - "m"''- (4.4) 

The extrinsic curvature tensor of the boundary is given by 

K'"' = /i^/i^V^n'^. (4.5) 

Its explicit expression in terms of the harmonic variables used here is given in App. 151 

The inhomogeneous boundary data q must be constrained if the Einstein equations are to be satisfied. The 
subsidiary system of equations which governs the constraints requires a slightly more complicated treatment than 
in the previous section. For its analysis, we represent the constraints by C = (C^,Ca = gafiC^). Since — 
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{C^rix — Can"-)^,^ + Cag"^^, it is clear that C = is equivalent to C = 0. It is straightforward to show, beginning 
with H2.11|l . that C satisfies a wave equation of the homogeneous form 

g^''9^9,C + P''a^C + QC = 0, (4.6) 

where the elements of the square matrices P'^ and Q depend algebraically on the evolution variables u. This system 
can be converted to symmetric hyperbolic first order form. Uniqueness of the solution to the IBVP then ensures that 
the generalized harmonic constraints are satisfied in the evolution domain provided the initial Cauchy data satisfy 
C = dtC^ = and that C satisfies a maximally dissipative homogeneous boundary condition. 



A. The boundary system 

We choose as the maximally dissipative homogeneous boundary condition for C 

Vg^n^C^le = C^\b = (4.7) 
^ /i^^S^aiB^g^a^Cals = 0. (4.8) 



Referring to ljA16|l . the resulting boundary flux associated with each component of C satisfies the dissipative condition 
J-^ = 0. These boundary conditions, along with the constraints on the initial Cauchy data, guarantee that solutions 
of the reduced system satisfy Einstein's equations. For simplicity in working out their consequence, we require the 
gauge source terms to satisfy 

n„f"|a = 0, (4.9) 

as in H3.11|l for the homogeneous case, but we make no assumption about the behavior of the constraint adjustment 
at the boundary. 

The boundary condition (|4.7|l . together with (|4.9() . then implies = = so that daj^"^ + dxj^^ — 0. Using 
(|4.1|) and (|4.2|) . this reduces to 

q-- = ~j--\Bdaq'\ (4.10) 

which provides the constrained Neumann boundary data 9^^. 

The boundary condition (|4.8() contains the term 9^7"^°, which must be eliminated by using the reduced equations. 
We proceed by using (|4.7|l to obtain, at the boundary, 

/i>^(V,C"' + V''a) = -h';^C^V,np + h';^n'^Vf3C, (4.11) 
= -K^pC" + h;£,,C, - h^Cf^W^n^ (4.12) 
= -2K^pC'' + h;CnC,, (4.13) 

where K^^'^ is the extrinsic curvature introduced in 14. 5() . In the a;''-coordinates of the boundary, this implies 

(V,C^ + W^Ca) - ^2^KabC' + g^^qf^d^Ca + g^'^Cbdaq". (4.14) 



Assuming that the reduced equations H2.9(l are satisfied, the boundary constraint (|4.8|l can then be re-expressed as 

KnpGi = -^g: = -K^bC" + l-^adaq' - -^a:. (4.15) 

V 5 ^ V 9 

which contains no second cc-derivatives of the metric. Here h^npG^ is the "boundary momentum" l|B9p so that this 
equation can be re-expressed as 

(utiK', - SlK) + KabC - ^Cbdaq' + -^A^) \b = 0. (4.16) 



This is a system of equations intrinsic to the boundary which provides the means for prescribing constraint-preserving 
values for the boundary data 
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B. Constraint-preserving boundary data 

The boundary system H4.16() can be recast as a symmetric hyperbolic system which determines constrained values of 
the Neumann data g"* in terms of the free (boundary gauge) data q°- and the boundary values of the other evolution 
variables U := [7°*, 7^^, i9a;7^")]B which cannot be freely specified. In doing so, there is considerable algebra in 
expressing (|4.16() in terms of the evolution variables and boundary data. 

First, the covariant derivative has to be explicitly worked out, as in ifgTUIl where we now set 5"'' = K"-'' - h''^K. 
From (|4.16(l . we then find (on the boundary) 

^b{^/^S'^') = P^ (4.17) 

where 



~h 



= -V^ih'^'S'"' - -h'^'S'^ddhc - V^K'^'a + ^V^^^h'^'Cddm" - ^-^h-'At- (4.18) 
Here 

For the adjustment (|2.14|) . 

f^ab^x ^ -aig''''C^dbq°'. (4.20) 
For the adjustment (|TTB|) . h^^Al = 0. For the adjustment (|TTB|) . 



XX 



For Fock's adjustment (|2.i:^|l . 

1 



^ab^x ^ ^^5^C'a^f_ (4 21) 



^abj^x ^ ^gXXfja^^g^g^ (4 22) 



From ljB7p . the individual terms constituting S are 



gab ^ XJL_i^h'"'dcq' + h!"'dcq'' -^h'^'^dcq") 



(9"' + + - l^q^dxl^' - qVdxj^'') , (4.23) 



2\/-hB g 

where it is assumed that is replaced by its constrained value H4.10|l . We write 



^/Zf^S"^ ^ ^i-q"'' + T"''), (4.24) 



where 



T''^ = ^/^ik'^'^dcq^ + h'"'dcq'' -2h"'''dcq'') 

- {^ + q''q^)q'''' + q''q^dxl''^ + q*'q^dxl'"' (4.25) 

does not contain Substitution of these pieces into (|4.17|) gives 

dbq"^ = (4.26) 

where 

V ^—[- q^^dbg'"'' + afc(g'='=T"'') - 2P" ) . (4.27) 
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This system of partial differential equations, with principle part dtq"'^, governs the boundary data q"''' in terms of the 
quantities U and q"". Here (which appears in P") must be expressed in terms of evolution variables according to 

C = ^{db-f"'' + a^7"^) - (4.28) 

We formulate a symmetric hyperbolic system by setting 

qab ^ Qab _ ^^ab^^^Qcd^ (4 29) 

where rjab is the auxiliary Minkowski metric associated with the coordinates of the boundary. Equation (|4.26(l then 
becomes 

dbQ""' - ^v'^'db ivcdQ^'') = (4.30) 

Next we introduce the "2+1" decomposition of the boundary — {t, x"^) — (t, y, z) and the quantities (j) ~ -^Q" and 
qA ^ QtA rpj^g^ (|On|l gives the symmetric hyperbolic system 

M + dsQ^ = V'-^dt[ScDQ^''] 

dtQ^ + S^^dB<t> = + i<5^^9B [<5ci3Q^^] - a^Q^^ (4.31) 

for the variables and Q^. Assuming that the boundary values of U are known, (|4.31l) determines the constrained 
boundary values of </> and Q-^ in terms of the free boundary data {q"", Q^^). 

Constraint-preserving boundary data q for the reduced system is determined from the free data (q", Q^^), with 
obtained from H4.10|) and q"^ obtained from H4.29|l and H4.31|l . This data q leads to the maximally dissipative 
homogeneous boundary data (|4.7|) and (|4.8|) for the subsidiary system (|4.6ll which governs the constraints. Conse- 
quently, any solution of the IB VP for the reduced equations with this boundary data is also a solution of the Einstein 
equations (assuming that the initial Cauchy data satisfies the constraints). 

In the case of the homogeneous boundary data q = treated in Sec. lIIII H4.10I) and (|4.16|) are automatically satisfied 
and the IBVP is well-posed. The appearance of the quantities U in H4.10|l and H4.16|l prevent application of Secchi's 
theorem to prove the well-posedness of the inhomogeneous IBVP. However, these theorems do imply a well-posed 
IBVP for boundary data ^q linearized off a solution with homogeneous data, which provides background values of 
U. From the analytical point of view, such linearized results often provide the first step in an iterative argument to 
show that the full system is well-posed. From a numerical point of view, which is our prime interest here, it suggests 
that the evolution code should be structured so that the values of U are updated before updating the values of q in 
an attempt to force U into such a background role. 



V. FINITE DIFFERENCE IMPLEMENTATION 



A first differential order formalism is useful for applying the theory of symmetric hyperbolic systems but in a 
numerical code it would introduce extra variables and the associated nonphysical constraints. For this reason, we 
base our code on the natural second order form of the quasilinear wave equations which comprise the reduced harmonic 
system (|2.9f) . They are finite differenced in the flux conservative form 

2V^E^d^ig''''df3j^'n~S^'' = (5.1) 

where S^'^ is comprised of (nonlinear) flrst derivative terms in H2.9(l that do not enter the principle part. 

Numerical evolution is implemented on a spatial grid {xi,yj, zk) = {Ih, Jh, Kh), < (/, J, K) < N, with uniform 
spacing h, on which a field f{t, x^) is represented by its grid values f[i^j^K]{t) = fit, xi, yj, zk)- The time integration is 
carried out by the method of lines using a 4th order Runge-Kutta method. We introduce the standard finite difference 
operators D^i and D±i according to the examples 

Doxfi.j — -^[fi+i,j~fi~i,j) 

D+xfi,,j = -^{fi+i,j - fi,j) 

D-yfij = -r{fi,.j~fi,j-i)\ 
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the translation operators T±i according to the example 

T±xfi,.i = fi±i,j (5.2) 
and the averaging operators A±i and Aoi , according to the examples 

A±xfi,j - i(r±, + i)/,,j 

Aoxfij = 2^"^+^ ^ T^x)fi,j- 

Standard centered differences Dqi are used to approximate the first derivative terms in (|5.1|l comprising S^'^, except 
at boundary points where one-sided derivatives are used when necessary. 

We will describe the finite difference techniques for the principle part of H5.1|) in terms of the scalar wave equation 

5„(g"^5^<i>) = 0. (5.3) 

The principle part of the linearization of H5.1|l gives rise to (|5.3|l for each component of 7'"'. The non- vanishing 
shift introduces mixed space-time derivatives dtdi which complicates the use of standard explicit algorithms for the 
wave equation. This problem has been addressed in |24l [34. 36, 37] . Here we base our work on evolution-boundary 
algorithms shown to be stable for a model 1-dimensional wave equation with shift 01 • The algorithms have the 
summation by parts (SEP) property that gives rise to an energy estimate for (|5.3I) provided the metric g^^ is positive 
definite, as is the case when the shift is subluminal (i.e. when the evolution direction t^da = dt is timelike). Alternative 
explicit finite difference algorithms are available which are stable for a superluminal shift ^3^ l37| | . Although such 
superluminal algorithms have importance in treating the region inside a black hole, they have awkward properties at 
a timelike boundary and are not employed in the test cases treated here. 

The algorithm we use is designed to obey the semi-discrete versions of two conservation laws of the continuum 
system H5.3f) . These govern the monopolc quantity 

Q = - / (5.4) 
Jv 

and the energy 

E ^\j^{-9''^l+g''^,^^,j)dV, (5.5) 

where dV — dxdydz. By assumption, the t — const Cauchy hypersurfaces are spacelike so that —17" > 0. We also 
assume in the following that g"^^ is positive definite (subluminal shift) so that E provides a norm. Note that in the 
gravitational case there are 10 quantities Q°'^ corresponding to H5.4|l . which have monopole, dipole or quadrupole 
transformation properties depending on the choice of indices. 

For periodic boundary conditions (or in the absence of a boundary), H5.3() implies strict monopole conservation 
Q^t = 0; and, when the coefhcicnts of the wave operator are frozen in time, i.e. when dtg"^ = 0, H5.3|l implies strict 
energy conservation E^t ~ 0. In the time-dependent, boundary-free case, 

E,t = \j^{9f^,c.<^,p)dV. (5.6) 
which readily provides an estimate of the form 

E^t < kE (5.7) 
for some k independent of the initial data for Thus the norm is bounded relative to its initial value at f = by 

E < Eoe'^K (5.8) 
More generally, these conservation laws include flux contributions from the boundary, 

Q.t = f N,g'"^,^dS, (5.9) 

JdV 
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and 

E^t^- I ^dS, (5.10) 

JdV 

where 

^ = -$,t5'«iV,$,„ (5.11) 

and where the surface area element dS and the unit outward normal Ni = are normalized with the Euclidean metric 
Sij . The monopole flux vanishes for the homogeneous Neumann boundary condition 

<?'"^.^.a = 0. (5.12) 

The energy flux is dissipative for boundary conditions such that J- > 0, e.g. a homogeneous Neumann boundary 
condition or a homogeneous Dirichlet boundary condition $.4 = or superpositions g^^Ni^^a + k^,t = with fc > 0. 

It is sufficient to confine the description of the finite difference techniques to the {x, y) plane and, for brevity, we 
present the evolution and boundary algorithms in terms of the 2-dimensional version of 1)5. 3|l . 



A. The evolution algorithm 



We now discuss the evolution algorithm in the 2-dimensional case with periodic boundary conditions /o,j = /at.j 
and //.o = fi,N- We define the semi-discrete versions of 1)5. 4|l and l|5.5|l as 



N 



^h' J2 (-5"'i>,t-/'^o.$). (5.13) 
(/,j)=i 



and 



N 



E = h^ £, (5.14) 



where 



+ -^{A+ygyy){D+y<i>f + \{A^ygyy){D^y<^f 

+ g''y{Do^^)Doy^. (5.15) 

The energy E provides a norm on the discretized system, i.e. E — implies ^ t — D±i^ — (provided g^^ is 
positive definite and the grid spacing is sufficiently small to justify the continuum inequalities resulting from positive- 
definiteness) . 

The simplest second order approximation to (|5.3|) which reduces in the 1-dimensional case to the SBP algorithm 
presented in is 

W := -at(g"9t$) - ^tig''Do^'^>) - Do^{g''^t<^>) - = (5.16) 

where 

Vl<^> = ^D+,(^{A_,g--)D_,<P^+^D_,(^{A+,g-^)D+,^ 

+ lD + y({A^ygyy)D^y^)+lD^y({A + ygyy)D + y^ 



Do. g^^Doy^ ]+Doy[ g'=yDo..<f . (5.17) 
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The semi-discrete conservation law Q t = for the case of periodic boundary conditions follows immediately from the 
flux conservative form of W. In order to establish the SBP property we consider the frozen coefficient case dtg"^ — 0. 
Then a straightforward calculation gives 



1 

r 

+ ^D+y(^iA^ygyy)iD^y<P)T^y<^>^t^ + ^ , ( ff^'^^ ) ( $) T+, $ , , 
+ (^<i>tT^,{g"Woy^) + g^y{Doy^)T^.^^t 

+ ^D+y (^<i>tr_^(5-^Z?o.$) + s'n^o.*)^-^*,*) . (5.18) 

Because each term in H5.18II is a total D±i, the semi-discrete conservation law E,t — follows (for periodic boundary 
conditions) when the evolution algorithm = is satisfied. When the coefficients of the wave operator are time 
dependent, an energy estimate can be established analogous to (|5.8|l for the continuum case. 

We also consider a modification of the algorithm H5.8|l by introducing extra averaging operators according to 

W := -dt{g''dt^) - dt (^g*'I?o.<f) - (^{A+,g''){A+A^)^ - DI<^ = 0, (5.19) 

with 



^D + y({A^ygyy)D^y^'] + ^D^y{ ( A + , ff^^ ) + , $ 



^ -D+4 {A^,,g--)D^,'S> ) + -D^4{A+,,g-^)D+^'P 



+ D-^i^{A+,g''y)iA+,Doy'S>)j+D^yi^{A+yg^y)iA+yDQ,'S>j. (5.20) 

It is easy to verify that W = W + 0{h?) and that both are constructed from the same stencil of points. Although W 
does not obey the exact SBP property with respect to the energy (|5.15|l . the experiments in [sj show that it leads to 
significantly better performance than W for gauge wave tests with periodic boundary conditions. 

For the time discretization, we apply the method of lines to the large system of ordinary differential equations 

= ^A*,i -f -1b*. (5.21) 



obtained from the spatial discretization. Introducing 

h 



*.t - tT, (5.22) 



we obtain the first order system 

T\ 1/BA\/T 



* ;^ h\l Q ) \ ^ 



(5.23) 

We solve this system numerically using a 4th order Runge-Kutta time integrator. 



B. The boundary algorithm 

Again it suffices to describe the algorithm in the 2-dimensional case. In the absence of periodic symmetry, we 
modify the definitions of the monopole quantity (|5.13|) and the energy (|5.14|) to include contributions from the cells 
at the edges = {(0, J), {N, J), (/, 0), (/, N)} and corners = {(0, 0), (A^, 0), (0, A^), (A^, A^)} of the rectangular grid: 

JV-l 

Q = h^ (-5"** - 5"^0^$) + Q-'^aes + Q-orner. (5.24) 

(/,J) = 1 edges corners 
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and 

E = h? ^ ^ £ + ^ ' E^dges + ^ ' Ecorners-i (5.25) 
(7,J) = 1 edges corners 

with £ given again by H5.15|l . We carry out the analysis of the semi-discrete conservation laws using the VF-algorithm 
(E3S|). 

We can isolate the contributions from the edges by retaining periodicity in the y-direction, so that the boundary 
consists of two circular edges (with no corners) at / = and I = N. We assign the contribution 

Qi=N = Y E - ff"** - 5*"^-x<i> - g'^Doy^ ) (5.26) 

,7=1 ^ 

to l|5.24(l from the I = N boundary, with the analogous contribution from / = 0. Assuming that W = a.t all interior 
points, a straightforward calculation yields 

Q,t = hY,{^Wo,j-No,j + ^WN,.j+MNj) (5.27) 
j=i 



where 



+ 5"^^oy$ + Ao4g^yDoy^)^ + y (^Doyig'^y D+^D^^^) + {dtg'^)D+^D-.,<i>^ . (5.28) 

Here, A/jv = is a second order accurate approximation to the Neumann condition H5.12|l . evaluated at the I — N 
boundary. Because Afpj involves ghost points outside the computational grid, the Neumann boundary condition must 
be implemented in the form 

'^Wo,j-Afo,j - 0, (5.29) 

which, by the above construction, involve only interior and boundary points. (This is why the term 
DQy{g'-^y D^xD-x^), which integrates to on the boundary, is included in H5.28|l .') After reducing (|5.29() to first 
order in time form as in H5.23|l . it provides the Neumann update algorithm for the boundary points when the interior 
is evolved using the W algorithm (|5.17|l . Correspondingly, when applying the W algorithm H5.19|l to the IB VP, we 
update the boundary using 

^Wn,.j+Mn,j = 
^Wq,.j-Nq,j = (5.30) 

where now 

^ ^ \ [iA+xgnD+x-^ + {A^^gno-.'^ + ^(.g*" + AW")r+x<i>t + i(ff*"' + ^_,g*^')r_,<i>t 

+ {A+xg-''')A+xDay^ + (A_,.g-«)A_,i?oy$ 

^ D^y[{A+yg^y)A+yD+xD-x^ + (9tg*")i?+.i?-.$ + i(i^+,I?_,g*")$t^ . (5.31) 
In the case of inhomogeneous Neumann boundary data g, H5.29|l takes the form 

^Wn,.j +Mn.j = g'^'-'qN 
^Wo,J-^fo,J = -g^^go. (5.32) 
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We assign the contribution to the energy from the I = N boundary to be 

,2 ^ / 1 1 

+ \{A+ygyy){D+y^f + ^{A_ygyy){D_y<i>r 

+ g"^p_,$)i^oy$), (5.33) 

with the analogous contribution from / = 0. Assuming that = at all interior points and that the coefhcients of 
the wave operator are time independent, this leads to the discrete energy- flux conservation law 

N 

E^t ^ -hY^i^oj + ^Nj) (5.34) 
J=i 

with 

^nj^-Im^W+U)) (5.35) 



'2 



N,J 



and 



^oj^ iM^W^M)] . (5.36) 



'2 



Any discrete boundary conditions such that J-q^j and J-'n,j are non-negative retain the SBP property of the W- 
algorithm. In particular, this holds for the Dirichlet or Neumann boundary conditions, for which = in the 
homogeneous case. We also consider the Sommerfeld-type algorithm (with Sommerfeld data q). 



V2 V -9tt J N, j 



Qn 



^W-M+.j^^t] ^ qo, (5.37) 



2 V -9tt 



0,J 



for which ^ is strictly positive in the case of homogeneous data. 

The algorithm for the corners follow from the same SBP approach. Since we restrict our applications in Sec. lVll to 
tests with smooth boundaries, we will not present the details. 

C. Dissipation 

It is simplest to add dissipation by using the Euclidean Laplacian with the centered approximation 

Vl = D+,D^, + D+yD^y + D+,D^,. (5.38) 
Then, for the 2-dimensional case with periodic boundary conditions, the modification 

W + eh^VlVl<i>^t (5.39) 
of the W evolution algorithm H5.16|l dissipates the energy H5.14|) according to 

N 

E^t^E^t + eh^ {'Dl'^.tf- (5.40) 

The dissipation H5.4()|l leaves undisturbed the semi-discrete energy estimate governing the stability of the W evo- 
lution algorithm. We present a similar SBP result for the boundary algorithm. It suffices to illustrate the approach 
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in the 1-dimensional case. It is the semi-discrete analogue of adding a dissipative term to the wave equation, in the 
half-plane < x, in the symbolic form 

5t(-.9"*.0 - + e{dt^,t + Six)dl<i>,t + d4Six)dl<i>,t)} (5.41) 

so that the Dirac delta function terms d{x) cancel the boundary contributions at x = to give 

^ dx'^>,tdt{-g"<i>^t)^ + (5.42) 
As a result, the energy dissipates according to 

poo 

dtE^dtE + e dx{dl^,tf. (5.43) 

In order to model this dissipation at the finite difference level, we use summation by parts techniques. A comprehen- 
sive discussion has been given in 23] . For the second order accurate approximation to the wave equation considered 
here, the following treatment suffices. 

On the interval Q < xi = Ih <1 — Nh we use the identity 

JV-2 

h J2 {fD-9 + {D+f)g}i = fN-i9N-2 - /2ffi (5.44) 

2 

to obtain 

hJ2{f{D+D^)'f}j = -hJ2{{D^f)D'_D+f}j + {DlD+fU^,fN^2-{D'-D+fhh 

2 2 

Af-2 

= hJ2 {{D+D^ff}i - {D^f )N-i{D+D^f)N-2 + (D^fUD+D^f)^ 

2 

+ {DlD+f)N-lfN-2 - {D^_D+f)2h 
N-2 

= hY,{{D+D^ff}i 

2 

r (fN-2D+D^fN-i - fN-iD+D^fN-2 + hD+D^h - fiD+D^f^ . (5.45) 



h ^ 

Thus we are led to the finite difference analogue of H5.41|) 

{dt{-9*'dt<S>)}i ~^dt{~g''dt<i>)j + eh^{D+D^fdt<S>i, 3<I<N-3 

at(-g"9t$), + e{h^iD+D^fdt<^>N-2 - hDldt<^>N-2} , I = N-2 
(9f(-g"9t$)^ + ehDldt<i>N-i} , / = iV - 1 
-.at(-g"9t$),, I^N, 1^0 

-> dt{-g''dt<S>)j + e{h\D+D^fdt<^2 - hOldt^^} , / = 2 
-^dt{-g''dt'i>), + ehDldt^i, 1 = 1. 

This dissipation implies (in the ID case where E = h'Y^E) 



(5.46) 



N-2 



E^t ^ E^t + eh^ J2 (D+D-dt^i f. (5.47) 

2 

Alternatively, we can extend the algorithm by one more point: 

{dt{-g'*dt<S>)}i '^dt{-g''dt^), + eh\D+D^fdt^i, 2<I<N-2 

-> dt{~g''dt<f)j + e{h\D+D^fdt<^N-i - hDldt<i>N-i} , I = N-1 
-^dti~g''dt<^>)j + ehD^_dt^N}, I = N 

dt{-g''dt<^>)j + e{h^{D+D^fdt<^>i - hDldt<^>i} , 1=1 
^ dti~g*'dt<S>)j + ehDldt^o, 1 = 0, (5.48) 
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with the dissipation now implying (in the ID case) 



(5.49) 



1 



Additional dissipation can be added at the boundary by using the identities 

KD^fN? = fND^fN - Jn-iD+Jn-i 



(5.50) 



and 



foD+fo + hD_h 



(5.51) 




dt{-g''dt^)j - eshD+dt^N-i , / = - 1 



dt{-g''dt^)j + eBhD^dt^N , I = N 
dt(-g''dt<^)j + eBhD^dt^i, 1 = 1 
dt{-g'*dt<l>)j~eBhD+dt'l>o, / - 0. 



(5.52) 



This results in the dissipative effect on the energy (in the ID case) 



(5.53) 



D. Implementing the constraint-preserving boundary system 



For a boundary located at a; = const, we have shown that the evolution-boundary algorithm is constraint-preserving 
if the boundary data (|4.2(l are provided and the remaining boundary values are updated as follows: 

1. The Dirichlet boundary data = j^a^^xx jg fj-ggj^y prescribed. 

2. The Neumann boundary data q^^ = g^(9p7^^ is determined by (|4.10() . 

3. The Neumann boundary data q""^ = q^df^-f'^'' is computed in terms of three free functions from the boundary 
system (|4.16ll . recast in the symmetric hyperbolic form H4.31|l . 

4. dtj^'^ is updated on the boundary using the constrained boundary data; then 7'*" is updated. 

The first two items involve minimal computation. The update algorithm ensures that the components 7^^ and 7°'' 
are known at the current time level. Then the Dirichlet data g° determines the remaining components 7°^ (by item 
1) and the Neumann data q^^ (by item 2). 

In describing the implementation of the third and fourth items, we consider for the present purposes a smooth 
boundary consisting of two toroidal faces at a; = and x = 1 (periodic boundary conditions in the y and z directions). 
Again, the grid points are labeled by < (/, J, K) < N. The updates of the boundary values 7^^ and terms 
of their Neumann boundary data is based upon the generalization of the scalar update scheme H5.32|l to include the 
non-principle-part terms of the gravitational equations. For example, the component of H5.32|l corresponding to 7°'' 
generalizes at the I = N boundary to 




= 0, 



(5.54) 





(5.55) 
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as in H5.28|l . 

The constrained Neumann data q""^ is obtained from the first order symmetric hyperbohc boundary system H4.31|l . 
which is solved by the method of hnes using centered spatial differences and a 4th order Runge-Kutta time integrator. 
The required data for this system (where — (y, z)) is related to the boundary data g"'' by (|4.29|) . or its inverse 



This gives 



Q-^_Qyv ^ q---qyy, 

Q^^ + Qyy = 2q'' ~ {q'"^ + qyy). (5.57) 
Integration of (|4.31|) then determines the constrained data 

QV = Qty ^ qty 

= Q*^ = (5.58) 

In the update of the Neumann components 7"'' via (|5.54|) . the inhomogeneous term g^^q°-^ introduces numerical 
error (from the evolution of g^^) which effects the exact nature of the semi-discrete multipole conservation laws 
satisfied by the principle part of the system. We avoid this source of error for the components corresponding to Q^^ 
by prescribing data for g^^Q^^ (rather than Q^^); similarly, we prescribe g^^q^^. In this way, the principle part of the 
system obeys exact inhomogeneous versions of the semi-discrete conservation laws for the components corresponding 
to Q^^. 



VI. TESTS OF THE EVOLUTION-BOUNDARY ALGORITHMS 



We compare the performance of various versions of the evolution-boundary algorithm using the AppleswithApples 
gauge wave metric 

ds"^ = (1 - H){-dt^ + dx^) -f dy^ + dz^, (6.1) 
and a shifted version given by the Kerr-Schild metric 

ds^ = -dt^ + dx^ + dy^ + dz^ + Hk^.kpdx^'dx^ , (6.2) 

where, in both cases, 

H = H{x-t)=A sin (^^^^ j , (6.3) 

and 



fca =<9„(x-i) = (-1,1,0,0). (6.4) 

These metrics describes sinusoidal gauge waves of amplitude A propagating along the x-axis. In order to test 2- 
dimensional features, we rotate the coordinates according to 

x=^ix'-y'), y = l=(x' + y'). (6.5) 

which produces a gauge wave propagating along the diagonal. 

The results of these gauge wave tests for periodic boundaries, i.e. a 3-torus without boundary, have been 
reported and discussed in [3l|. In the periodic case, it was found that the ID tests were very discriminating in 
revealing problems. The 2D tests were essential for designing algorithms for handling the mixed second derivatives 
in the Laplacian operator but they introduced no new instabilities of an analytic origin. As a result, in the ID tests, 
numerical error is channeled more effectively into exciting nonlinear instabilities of the analytic problem. 
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Here we open up the a;- axis of the 3-torus to form a manifold with boundaries at x = ±.5. Most of our 
tests are run, in both axis-ahgned and diagonal form, with amplitude A — 0.5. We have found that the smaller 
amplitudes, A — .01 and A ~ .1, of the original AwA specifications are not as efficient for revealing problems. Larger 
amplitudes can trigger gauge pathologies more quickly, e.g 5*' > (breakdown of the spacelike nature of the Cauchy 
hypersurfaces). Otherwise, the tests are carried out with the original AwA specifications. We choose wavelength d = 1 
in the axis-aligned simulation and wavelength d' = \/2 in the diagonal simulation and make the following choices for 
the computational grid: 

• Simulation domain: 

ID: a; G [-0.5, +0.5], y = 0, z = 0, d=l 

diagonal: x e [-0.5, +0.5], y E [-0.5, +0.5], 2 = 0, d' = V2 

• Grid: x„ — —0.5 + ndx, n = 0, 1 . . . 50p, dx — dy = dz = l/(50p), p = 1, 2, 4 

• Time step: dt = dx/i = 0.005//9 . 

The grids have N = 50p = (50, 100, 200) zones. (At least 50 zones are required to lead to reasonable simulations for 
more than 10 crossing times.) The ID tests are carried out for t — 1000 crossing times, i.e. 2 x lO^p time steps (or 
until the code crashes) and the 2D tests for 100 crossing times. 

For the case without periodic boundaries in the cc-direction, the boundary data are provided at a; = ±.5. For 
example, in the ID simulation, for a formulation with a Sommerfeld boundary condition on the metric component gtt 
the correct boundary data are 

{dt - d.,)gttU=^.5 = 2dtHi-.5-t) (6.6) 
[dt + d.,)gtt\cc=.5 = 0. (6.7) 

With this inhomogeneous Sommerfeld data, the wave enters through the boundary at x = —.5, propagates across 
the grid and exits through the boundary at x — .5. Inhomogeneous Dirichlet and Neumann data are supplied in 
the analogous way. Analytic boundary data are prescribed for all 10 metrics components except for the constraint 
preserving boundary algorithm, where only 6 pieces of analytic data are provided. 

The test results reported here are for the W algorithm (|5.19|l . which performed significantly better than the W 
algorithm in tests in the boundary-free case [sj. Also, harmonic gauge forcing terms were not found to be effective 
in the boundary-free gauge wave tests and we have not not included them in the present tests. (Gauge forcing is 
important in spacetimes where harmonic coordinates are pathological, e.g. the standard t-coordinate in Schwarzschild 
spacetime is harmonic but singular at the horizon.) 

We use the £00 norm to measure the error 

£ = ||$p-$,„a||oo (6.8) 

in a grid function <f>p with known analytic value $ana- We measure the convergence rate at time t 

Kt)=log,( ||^^'^°"°|| ), (6.9) 

11^4 ^ana 1 1 

using the p = 2 and p — A grids (A^ = 100 and — 200). It is also convenient to graph the rescaled error 

Sp=^\\^p~^ana\\^, (6.10) 

which is normalized to the p = 4 grid. 



A. Gauge wave boundary tests 

Simulation of the AwA gauge wave without shift H6.1|l is complicated by the related metric ^3] 

dsl = e^\l - H){-dt^ + dx"^) + dy"^ + dz^. (6.11) 

For any value of A, this is a flat metric which obeys the harmonic coordinate conditions and represents a harmonic 
gauge instability of Minkowski space with periodic boundary conditions. Since H6.11|l satisfies the harmonic conditions. 
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constraint adjustments are ineffective in controlling this instability. Long term harmonic evolutions of the AwA gauge 
wave with periodic boundary conditions were achieved by suppressing this unstable mode by the use of semi-discrete 
conservation laws for the principle part of the system. 

These semi-discrete conservation laws have been incorporated in Sec. in the harmonic evolution-boundary algo- 
rithm with a general dissipative boundary condition. Consequently, good long term performance should continue to 
be expected. This is especially true for Dirichlet or Sommerfeld boundary conditions, which do not allow the unstable 
mode (|6.11|l at the continuum level. However, Neumann boundary conditions allow this mode so that it poses a 
potential problem for the constraint preserving boundary algorithm described in Sec. IV Dl which combines Dirichlet 
and Neumann boundary conditions. Nevertheless, our expectations of good long term performance are borne out in 
all our tests of the gauge wave without shift. No explicit dissipation was added to the Ty-algorithm. 

Figures m 121 a-nd 121 plot the rescaled error Ep{t) in g^x for tests with amplitudes A — .5 and either 10 Dirichlet, 10 
Neumann or 10 Sommerfeld boundary conditions. The plots for the Dirichlet and Neumann boundary conditions show 
accuracy comparable to the corresponding tests with periodic boundary conditions reported in |3l[ . The convergence 
rates measured at i = 50 are rioD(50) = 1.876 for the Dirichlet case and rioAr(50) = 1.804 for the Neumann case. 
The Dirichlet and Neumann boundary algorithms supply the correct inhomogeneous boundary data for the gauge 
wave signal to leave the grid but the numerical error is reflected at the boundaries and accumulates in the grid. The 
Sommerfeld boundary condition allows this numerical error to leave the grid and gives excellent long term performance, 
as evidenced by the clear second order convergence displayed in Fig. O throughout the 1000 crossing time run. The 
convergence rate rios(50) — 2.000 was measured at t = 50. 
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FIG. 1: Plot of the rescaled error £p{t) in Qxx for the ID gauge wave with 10 Dirichlet boundary conditions using the bare 
W algorithm (no constraint adjustment or dissipation) obtained with — 50, 100 and 200 grid zones. The horizontal axis 
measures crossing time t. The instabilities are kept under control during the entire 1000 crossing times run. 
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FIG. 2: Plot of the rescaled error £p{t) in g^x for the ID gauge wave as in Fig. but with 10 Neumann boundary conditions. 

Test results using 3 Dirichlet boundary conditions (for the components 7^") and seven Neumann boundary con- 
ditions (for the components 7^^ and 7"''), with boundary data supplied by the analytical solution, were practically 
identical to results for the 10 Neumann boundary conditions. When these seven pieces of Neumann boundary data 
was generated by the constraint preserving boundary system, as detailed in Sec. lVDl practically identical results were 
again found. This is because no appreciable constraint violation is excited in the shift-free gauge wave test. For the 
same reason, constraint adjustments have no significant effect. 
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FIG. 3: Plot of the rescaled error £p{t) in g^x for the ID gauge wave as in Fig.Qbut with 10 Sommerfeld boundary conditions. 
Excellent performance with clean second order convergence is manifest for 1000 crossing times. 



Two-dimensional tests of the diagonally propagating gauge wave were also in line with expectations. The most 
important background information for the design of physically relevant boundary algorithms using CCM comes from 
the case of 10 Sommerfeld boundary conditions and the case of boundary conditions supplied by the constraint 
preserving (CP) boundary algorithm. In 2D tests of these algorithms, we found the respective convergence rates 

rios(lO) = 2.0005 , rios(lOO) = 2.0004 
rcp(lO) = 2.0482 , rcp(lOO) = 1.0338 

measured at 10 and 100 crossing times. The Sommerfeld case maintains clean second order convergence. At 10 
crossing times the constraint preserving system also shows clean second order convergence but at 100 crossing times 
the convergence rate drops to first order due to accumulation of error. The constraint boundary system (|4.31|1 involves 
two derivatives of numerically evolved quantities in order to provide a complete set of Neumann data to the code. 
This introduces high frequency error in the Neumann data. We have not experimented with ways to dissipate this 
source of error in the boundary system. 



B. Shifted gauge wave boundary tests 
Simulation of the shifted gauge wave 1)6. 2|l is complicated by the related exponentially growing Kerr-Schild metric [s^ 

dsl = -dt^ + dx^ + dy^ +dz^ + (^H -1 + e^*^ kakpdx°'dx'^ , (6.12) 

where 

Ad /2Tr(x-t)\ 



Air 

Although this metric does not solve Einstein's equations, for any value of A it satisfies the standard harmonic form 
(|2.4() of the reduced Einstein equations, i.e. the equations upon which the numerical evolution is based. Numerical 
error in the shifted gauge wave test should be expected to excite this instability. As a result, in shifted gauge wave 
tests with periodic boundary conditions ^i|j the constraint adjustment (|2.14|l or (|2.15|) was necessary to modify the 
form of the reduced system in order to eliminate the solution H0.12|l . 

Figures 01 |5l and plot the rescaled error £p{t) in g^x for shifted gauge wave tests with amplitude A — .5 using the 
bare W algorithm (no constraint adjustment or dissipation). The results in Fig . IH for 10 Dirichlet boundary conditions 
are much less accurate than test results for the bare algorithm reported in |31j for periodic boundary conditions. The 
periodic tests run about 4 times longer with the comparable error. The results for 10 Neumann boundary conditions 
shown in Fig. are poorer yet by an additional factor of 4 when compared to the periodic tests. These results can be 
attributed to the periodic motion of the boundary introduced by the shift. As explained in more detail in p^. the 
numerical noise can be blue shifted as it is reflected at the moving boundaries, leading to rapid build up of error. The 
results show second order convergence only at early times. For the 10-Dirichlet test we measured the convergence rate 
''iod(IO) = 2.077 at 10 crossing times. For the 10-Neumann test, we measured the convergence rates riojv(l) = 2.013, 
''ioAr(5) = 2.760 and rioAr(lO) = 3.978, at 1, 5 and 10 crossing times. The high convergence rate at 10 crossing 
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times is misleading since instabilities have already dominated the N = 100 run (the coarser grid in the convergence 
calculation). 

Figure El shows that 10 Sommerfeld boundary conditions gives very good long term performance for the bare 
algorithm. The accuracy after 1000 crossing times is much better than found in *3l] for the bare algorithm with 
periodic boundary conditions, as expected since the Sommerfeld conditions allow numerical error to leave the grid. 
We measured the convergence rate 

rios(50) = 2.413 (6.14) 

at 50 crossing times. 




FIG. 4: Log plot of the rescaled error £p{t) in g^x for the ID gauge wave with shift for the bare W algorithm with 10 Dirichlet 
boundary conditions. Although clean second order convergence is measured at early times, the simulation goes unstable in less 
than 100 crossing times. 




FIG. 5: Log plot of the rescaled error £p{t) in g^x for the ID gauge wave with shift as in Fig.|3]but with 10 Neumann boundary 
conditions. Again the runs converge at early times but now the unstable mode is evident at less than 30 crossing times. 

As we have already pointed out for the case of periodic boundary conditions, although clear second order convergence 
was established for short run times, good long term performance for the shifted gauge wave tests was only possible 
when constraint adjustments were used to control the unstable mode (|6.I2|I . The situation is similar in the presence of 
boundaries. Figure shows the dramatic improvement in performance obtained by using the adjustment H2.15|l . with 
02 = 1, for ID tests with 10 Dirichlet boundary conditions, with no dissipation added. The figure shows good second 
order convergence, with only a small, non-growing error at 1000 crossing times. The convergence rate ?'io£)(50) = 2.018 
was measured at 50 crossing times. For the same adjustment, with no dissipation. Figure |S1 shows the test results for 
10 Neumann boundary conditions. Again there is very good accuracy for 1000 crossing times, although the error is 
larger than the Dirichlet case because the boundary conditions now allow the instability (|6.12() to be excited. The 
convergence rate rioAf(50) = 2.877 was measured at 50 crossing times. 

The ID tests with 10 Sommerfeld boundary conditions already gave good results for the bare algorithm which 
constraint adjustment or dissipation do not substantially improve. Figure El shows the rescaled error for 2D shifted 
gauge wave tests with the 10-Sommerfeld boundary algorithm. The convergence rates 



rios,2D(10) = 2.096 , rios,2D(100) = 2.680 



(6.15) 
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FIG. 6: Plot of the rescaled error £p{t) in g^^ for the ID gauge wave with shift as in Fig. 2]but with 10 Sommerfeld boundary 
conditions. The error has only a slow linear growth in time. 



were measured at 10 and 100 crossing times. 
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FIG. 7: Plot of the rescaled error £p{t) in g^x for the ID shifted gauge wave test with 10 Dirichlet boundary conditions as in 
Fig. |1] but when the adjustment H2.15|l is used. The dramatic improvement compared with Fig. |l]is evident. There is no long 
term error growth during the 1000 crossing time run. 
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FIG. 8: Plot of the rescaled error £p{t) in g^x for the ID shifted gauge wave test with 10 Neumann boundary conditions as in 
Fig. |S] but when the adjustment 12.1511 is used. Again, the adjustment leads to dramatic improvement. 



The constraint preserving boundary algorithm described in Sec. IIVI requires 3 Dirichlet boundary conditions (for 
7^°) and 7 Neumann boundary conditions (for 7^^ and 7"''). As shown in Fig llOl when the 3 Dirichlet and 7 
Neumann pieces of boundary data are suppUed by the analytic solution the results for the shifted gauge wave test 
are very poor due to the early excitation of an unstable error mode. The rescaled error plotted in Fig. 1101 shows that 
good performance is maintained for less than 8 crossing times. The constraint adjustments H2.14|l - H2.1(j|l . as well 
as other adjustments considered in [sj, do not lead to improvement. In particular, the adjustment (|2.15|) . which 
works remarkably well in suppressing instabilities in the 10-Dirichlet and 10-Neumann algorithms, fails to be effective 
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FIG. 9: Plot of the rescaled error Sp{t) in Qxx with 10 Sommerfeld boundary conditions as in Fig.inibut now for the 2D shifted 
gauge wave test. As in the ID test, the error has only a slow linear growth in time. 
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FIG. 10: Log plot of the rescaled error £p{t) in gxx for the ID shifted gauge wave test with 7 Neumann and 3 Dirichlet 
boundary conditions and amplitude A = 0.5. The instabilities which appear in less than 8 crossing times cannot be controlled 
by dissipation or by the constraint adjustments considered here. 

when these algorithms are mixed. This seems related to the fact that a Kerr-Schild pulse reflected by 10-Dirichlet 
or 10-Neumann conditions results in a Kerr-Schild pulse traveling in the opposite direction. When the Dirichlet and 
Neumann conditions are mixed, the reflected pulse has opposite signs for the Neumann and Dirichlet components. As 
a result, the reflected pulse no longer has the necessary Kerr-Schild properties for which the adjustment H2.15|l was 
designed and a new unstable mode is excited. The effect arises from the nonlinear coupling between the Dirichlet and 
Neumann components and does not arise in the AwA gauge wave test where the Dirichlet components vanish. The 
new unstable mode has long wavelength so it cannot be controlled by dissipation and it satisfles Q ~ (see (|5.4|) ') so 
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FIG. 11: Log plot of the rescaled error £p{t) in g^x for the ID shifted gauge wave test with 7 Neumann and 3 Dirichlet boundary 
conditions as in Fig. llOl but with amplitude A = .1. 
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FIG. 12: Log plot ol the rescaled error £p(t) in g^x for the ID shifted gauge wave test with 7 Neumann and 3 Dirichlet boundary 
conditions as in Fig. llfll but with amplitude A — .01. The rapid attenuation of the growth rate of the instabihty with amplitude 
is evident from comparison with Fig's ITUl and im 

that it cannot be controlled by the semi-discrete conservation law Q,t = 0. 

The effect of the constraint-preserving boundary system on the shifted gauge wave test depends upon the perfor- 
mance of the underlying 3-Dirichlet, 7-Neumann boundary algorithm so it also does not lead to good performance for 
amplitude A — .5. For smaller amplitudes, the instability from the Dirichlet- Neumann mixing is weaker. Figures ITTI 
and 1121 show how the rescaled error behaves as the amplitude is lowered to ^ = .1 and A — .01, when the Dirichlet- 
Neumann data is supplied analytically. For the N — 200 grid and amplitude A — .01, reasonable performance is 
maintained for 500 crossing times. 

Only slight improvement is obtained when the constraint preserving boundary system is applied to the 3-Dirichlet, 7- 
Neumann algorithm. This is apparently because the troublesome instability is a constraint preserving mode. Figures 
1131 and 1141 compare the £cx) error norms for constraint violation with and without the application the constraint 
preserving boundary system for the shifted gauge wave with amplitude A = A. The figures show that both the 
and C* components of the harmonic constraints are satisfied to much higher accuracy when the constraint preserving 
boundary system is used. However, the unstable mode still grows on the same time scale and eventually dominates 
the simulation. 




FIG. 13: Log plots of the error £p{t) in the constraint for the ID shifted gauge wave test with amplitude A — A with 7 
Neumann and 3 Dirichlet boundary conditions on an = 200 grid. The curve (• ■ •) is the error when the boundary data is 
supplied analytically and the curve ( ) is the error when the constraint preserving boundary system is applied. 



VII. SUMMARY AND IMPLICATIONS FOR CAUCHY-CHARACTERISTIC MATCHING 

We have formulated several maximally dissipative boundary algorithms for well-posed and constraint preserving 
versions of the general harmonic IBVP. The algorithms incorporate SEP and other semi-discrete conservation laws 
for the principle part of the system. In boundary tests based upon the AwA gauge wave, we have demonstrated that 
these techniques give rise to excellent long term (1000 crossing times), nonlinear (amplitude A = .5) performance for 
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FIG. 14: Log plots of the constraint error £p{t) as in Fig. 1131 but for the component C'. Although the constraint preserving 
boundary system ( ) leads to much less constraint violation, the time scale of the unstable mode is unaffected. 

all boundary algorithms considered. This is an especially challenging test. A high amplitude periodic wave enters 
one boundary, crosses the computational domain and exits the other boundary. The numerical error modes, which 
are continously being excited, are trapped between the boundaries by the reflecting Neumann or Dirichlet conditions. 
However, this causes no serious problem because the chief unstable mode is suppressed by the conservation laws and 
the growth of error is limited by the dissipation inherent in the system. 

The shifted gauge wave boundary test proved even more challenging, as expected from the blue shifting of the 
error resulting from reflection off the oscillating boundaries. However, 10 Sommerfeld boundary conditions (on the 
components of the metric) continued to give excellent long term, nonlinear performance. Furthermore, with the use of 
the constraint adjustment H2.15|l . 10 Dirichlet or 10 Neumann boundary conditions also continued to give good long 
term, nonlinear performance. However, nonlinear instabilities, which did not respond to constraint adjustment, were 
excited by the mixture of 3 Dirichlet and 7 Neumann boundary conditions. This led to much poorer results and rea- 
sonable long term performance required amplitudes A < .01. When the 3-Dirichlet, 7-Neumann boundary algorithm 
was amended to include the constraint preserving boundary system, the performance was only slightly improved. 
The results for this test problem indicate that the mixing of Dirichlet and Neumann boundary conditions excites a 
nonlinear instability of the analytic problem which cannot be controlled by the numerical techniques considered here. 

The above results focused on the mathematical and computational issues underlying accurate long term simulation of 
the IBVP but they also have bearing on the important physical issue of providing proper boundary data for waveform 
extraction by means of CCM. In the absence of analytic boundary data, as supplied in the gauge wave testbeds, 
constraint preserving boundary data must be obtained either by using a constraint preserving boundary system or by 
matching to an exterior solution which accurately satisfies the constraints. The inaccuracies which we have found in 
applying the constraint preserving boundary system would only be tolerable in matching with a Cauchy boundary in 
the weak field regime, i.e. far from the sources. On the other hand, the use of 10 Sommerfeld boundary conditions 
shows robust long term accuracy even in the strong field regime and would allow the economy of a boundary close to 
the sources. This is a strong recommendation for a matching algorithm based upon Sommerfeld boundary conditions, 
with constraint preserving Sommerfeld data supplied by matching to an exterior characteristic solution. 

Constraint preservation with a Sommerfeld boundary algorithm in CCM is greatly facilitated by using harmonic 
evolution. In CCM, there is a Cauchy evolution region with outer boundary Bq and a characteristic evolution 
region, extending to I+, with inner boundary Bi inside the Cauchy boundary Bq- A global solution is obtained by 
injecting the necessary Cauchy boundary data on Bq from the characteristic evolution. In doing so there are two 
concerns: (i) injection of the Cauchy data in the correct gauge and (ii) injection of constraint-preserving Cauchy 
data. Harmonic evolution simplifies dealing with both of these concerns because the gauge conditions and constraints 
reduce to wave equations for the harmonic coordinates x". By extracting data for on the inner characteristic 
boundary Bi, the harmonic coordinates can be accurately propagated by characteristic evolution to the injection 
worldtube Bq as solutions of the wave equation. This supplies the necessary Jacobian between the characteristic and 
Cauchy coordinates to ensure that the above two concerns are properly handled. Furthermore, since such data is 
constraint-preserving (up to numerical error), it is possible to inject the data in (inhomogeneous) Sommerfeld form. 
As evidenced by the test results in Sec. IVII this offers a very robust way of injecting boundary data which obeys all 
the constraints of the harmonic system, while avoiding the computational error introduced by solving the boundary 
constraint system. This was part of the strategy employed in 6] in successfully implementing long term stable CCM 
in linearized gravitational theory. The results of this paper suggest that this is also a promising approach for carrying 
out CCM in the nonlinear case. 



28 



APPENDIX A: IBVP FOR A SCALAR FIELD IN A CURVED BACKGROUND SPACETIME 

The well-posedness of the IBVP for the scalar wave equation with shift can be estabhshed by standard techniques 
by reducing the evolution system to symmetric hyperbolic form. We consider the principle part of the scalar wave 
equation written in the form 

gt^'d^d.'P = 0. (Al) 

We reduce this to the first order symmetric hyperbolic system A^dtu + A^diU — S, following a treatment in by 
introducing the auxiliary variables T = dt^, X = dx^, y = 9y<I>, and Z = 9^$. Then in terms of the 5 dimensional 
column vector u — ^{^, T, X, y, Z), the matrices A" are given by 

1 

^* = I -5** I (A2) 



g- 



jk 



and 





A' = \ ~2g^' ~gi^ | (A3) 
-g'^ 

and S — '^{T, 0, 0, 0, 0). In this first order form, the Cauchy data consist of uq = u\t=Q subject to the constraints 

{Cx,Cy,c,) := (x, y, z) - dy'^, = o. (A4) 

The evolution system implies that the constraints propagate according to 

dt{Cx,Cy,C,) = Q. (A5) 

The well-posedness of the Cauchy problem follows from the established properties of symmetric hyperbolic systems. 

We now consider the IBVP in the domain t > Q with a; < 0. The constraints propagate up the timelikc boundary 
at a; = and present no complication. The "boundary matrix" 

/O \ 

A"" ^ \ Q -25*^ -gJ^ (A6) 

yo I 

has a 3-dimensional kernel, with a basis corresponding to the vectors 

^(1,0,0,0,0), ^(0,0,-5^^3^^0), and ^(0, 0, Og^^). (A7) 
In addition, there is one positive eigenvalue and one negative eigenvalue 

A± = ±A - g"*, (AS) 

where 



X = ^{g^t)2 ^S^^gX^gXJ_ (A9) 

Thus precisely one boundary condition is required. The corresponding normalized eigenvectors are 

u± = I ^(0,TA + ff^^5^^^7^^5"^). (AIO) 

\/±2. A A± 

The homogeneous boundary condition Mu = may be cast in the form that u lies in the linear subspace (u+ + 
Hu- + Uq), where uq lies in the kernel. In this subspace, the flux 

T"" = - ^{u+ + Hu^ + Mo)A"^(u+ + Hu^ + uo) = A+ + A_ff^ (All) 
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satisfies the dissipative condition > provided 

< -A+ / A_ . (A12) 
The homogeneous boundary condition corresponding to a choice of H is 

Mu:=^(iJM+-M_)u = 0. (A13) 



The limiting case H = ^ — A+ / A_ leads to the homogeneous Dirichlet boundary condition 

= (A14) 



and the limiting case H = — ^ — A+ / A- leads to the homogeneous Neumann boundary condition 

= 0. (A15) 

As an alternative to this standard eigenvalue analysis of the maximally dissipative boundary condition, the simplicity 
of the scalar wave case allows a direct geometrical approach. A straightforward calculation gives 

J^^ = i ^uA^u = -(9t$).g^"9„$, (A16) 

which is identical to the energy flux obtained from the stress-energy tensor of a massless scalar field. It immediately 
follows that = for the above specifications of homogeneous Dirichlet and Neumann boundary conditions. More 
generally, J-^ > for boundary conditions of the form {g^^d^^ + Pdt^), with P > 0. An important case is the choice 



P=\-— (A17) 
9tt 



which leads to the homogeneous Sommerfeld-like boundary condition 



L^'a^* = g^'^a^* + J-^dt^ = 0, (A18) 
V 9tt 

where lies in the outgoing null direction in the (i", V"x) plane. 

The well-posedness of the IBVP for the scalar wave equation with any of the above choices of maximally dissipative 
boundary condition extends, by Secchi's theorem |4lj |. to the quasilinear reduced Einstein equations (|5.1() . The 
extension of the homogeneous boundary condition Mu = to include free boundary data q takes the form M{u—q) — 0. 
For example, the Neumann boundary condition ljA15|l takes the inhomogeneous form g"^da^ = q, where the boundary 
data q can be freely assigned. The well-posedness of the IBVP for the wave equation with inhomogeneous boundary 
data follows from the well-posedness of the homogeneous case by a standard argument. 

APPENDIX B: EXTRINSIC CURVATURE 

The extrinsic curvature tensor of the boundary at a; = (of the domain a; < 0) is given by 

j^^u ^ h'^h'^pV^nl^ (Bl) 

where 

K=S'i.~ n^n'^ (B2) 

is the projection tensor into the tangent space of the boundary and Ua — {1/ \/g^)'^ aX is the unit outward normal. 
In the 3+1 decomposition — {x°',x), hat is the intrinsic metric of the boundary, — 0, h^h = g^b and h^x — 
gxx — (l/ff^^^-) Also note that = det{hab) — g^^g- 

We set = ^g^^q^'- — \/g^{q°' , 1), where q° — j^^^ . In terms of the metric, 

K^^ = y^h'^h'-^^g-Pd^qP + gf^Pg^q- _ qPd,g»^). (B3) 
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In expressing this in terms of 7^*^, the identity 

dcg = -V^g^.dc.i'^'' (B4) 

is useful. For example, 

d.g^-" = -^(5.7'"^ - ^g'^-'gcpd.^r^). (B5) 

V 5 ^ 



As a result, 



and 



- -^(M^/j-qpa^7"/3 + /,M^„^^^qPa^^o/5)^ 

The boundary version of the usual momentum constraint implies 

Kn^C" = D^K'"' - h^^'^K), (B8) 

where Da is the covariant derivative associated with the boundary metric hab- We express this "boundary momentum" 
in the 3-dimensional form of the x"" — (t,x,y) coordinates: 

Kn^G"" = DbiK"'' - h^-^K) = 0. (B9) 

The identity 

^/^^DbSt = dbi^/^S'j - ^^^S'^dahc, (BIO) 
valid for any symmetric tensor S^'^, is useful in calculating the covariant derivatives entering the constraint. 
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